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Time-resolved electron transport in nano-devices is described by means of a time-nonlocal quan- 
tum master equation for the reduced density operator. Our formulation allows for arbitrary time 
dependences of any device or contact parameter. The quantum master equation and the related 
expression for the electron current through the device are derived in fourth order of the coupling to 
the contacts. It is shown that a consistent sum up to infinite orders induces level broadening in the 
device. To facilitate a numerical propagation of the equations we propose to use auxiliary density 
operators. An expansion of the Fermi function in terms of a sum of simple poles leads to a set of 
equations of motion, which can be solved by standard methods. We demonstrate the viability of 
the proposed propagation scheme and consider electron transport through a double quantum dot. 
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I. INTRODUCTION 

One of the most successful techniques for investigating 
artificial nano-structures, such as integrated molecules or 
gated quantum dots, is the measurement of the electron 
flow through the structure from a source to a drain con- 
tact. In this field, there has been an enormous progress in 
the experimental realization since the first measurements 
of electric currents through molecules [l|, H| or quantum 
dots [1, H| more then a decade ago. Thereby, it has 
been successfully demonstrated that current measure- 
ments provide information equivalent to the traditional 
spectroscopy of atoms or molecules in the gas phase |5j. 
In recent years, the possibility to inject ultra-fast volt- 
age pulses into nano-devices paved the way to perform 
transport measurements beyond steady-state paradigms. 
This approach is particularly promising for two reasons: 
Firstly, transient relaxation processes can be measured 
directly in the time domain Secondly, the ex- 

ternal control of the pulses can be used to manipulate 
the electronic dynamics in quantum-dot systems |8l4l3lj . 
which may eventually be used to realize quantum gates 
14l | . In either case an explicit time-dependent theo- 
retical description of the ultra-fast dynamical response 
is mandatory for interpreting time-resolved data. Fur- 
ther, it would possibly allow to develop new manipula- 
tion schemes for nanoscale devices. 

Addressing these challenges there has been consider- 
able progress towards having real-time methods for the 
description of electron transport. Many of these are 
based on non-equilibrium Green functions (NEGF) [l5r - 
22j ■ Hereby, electron-electron interaction is either fully 
neglected or approximately taken into account by means 
of self-energies obtained from many-body perturbation 
theory [2l| . An alternative approach for strongly inter- 
acting systems is given by considering the time evolution 
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of the (reduced) many-particle density matrix [23W30I ]. 
As in the theory of open quantum systems [3l[ , one stud- 
ies the reduced density matrix describing exclusively the 
degrees of freedom of a particular sub-system (here the 
nano-device) which is embedded in a larger system con- 
taining an environment (here the contacts). Instead of 
the bosonic environment of a system coupled to a heat 
bath, the nano-device is coupled to thermal fermionic 
reservoirs through source and drain contacts. The tun- 
nel coupling of these contacts to the system leads to an 
exchange of electrons with the reservoirs, and thus to an 
electric current through the device. This current, as well 
as related quantities like level occupations, can be deter- 
mined by means of the reduced density matrix [27| . Its 
time evolution is given by so-called (generalized) quan- 
tum master equations (QME), in which typically the cou- 
pling to the environment is incorporated in a perturba- 
tive description, often in second order. 

In this article, we present a propagation scheme which 
allows for an efficient real-time description of nano- 
systems with arbitrary time-dependences of any parame- 
ter of the system and the contacts, such as level energies, 
tunnel couplings, or voltages. We discuss a fourth-order 
method, which serves two purposes: Firstly, it shows cor- 
rections beyond the "standard description" with second- 
order schemes [H, [27], . Secondly, it devices a method 
for summing up certain terms up to infinitely many or- 
ders. Thereby, the broadening of system levels due to the 
coupling to the contacts is properly taken into account, 
which is shown to be crucial for obtaining meaningful 
results from solving QMEs in the time domain. 

Additionally, we propose the application of auxiliary 
density operators by which the calculation of the time 
evolution is reduced to the propagation of coupled dif- 
ferential equations containing only quantities with one 
time argument (27l . |32| . By means of a recently pro- 
posed partial fraction decomposition of the Fermi func- 
tion [33| the number of equations to be propagated is 
considerably smaller than in previous implementations 
of the auxiliary-operator method [27l . l28j . 
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In the following we define the setup of system, contacts 
and their coupling by giving the corresponding Hamilton 
operators. In Sect. Ql] we summarize the derivation of 
the QME for the reduced density matrix up to fourth or- 
der. We show how the current can be calculated in terms 
of reduced (i. e. system) quantities. Auxiliary operators 
are introduced in Sect. Mil We derive their equations of 
motion and give an expression for calculating of the elec- 
tric current. In Sect. II Vl we discuss the extrapolation to 
infinite orders. Finally, we present in Sect. [V] a numer- 
ical application of our propagation scheme to modeling 
double quantum dots. 



A. Setup 

We consider the usual threefold setup for the descrip- 
tion of tunneling problems, i. e. a device which is coupled 
to two electron reservoirs via tunneling barriers [3J, |35| . 
Accordingly, the total Hamiltonian can be split into three 
parts, 



H = H s + Hr 



H 



SR • 



The device Hamiltonian 



H, = 



£l(t)c\ci + ^ V lm(t)cjc m + Hi 



(i) 



(2) 



is characterized by discrete levels si which may be cou- 
pled through Vi m as, e.g., levels from different parts in 
a quantum dot system. In Eq. ([2j and below the opera- 
tors cj and ci denote the creation and annihilation of an 
electron in state Electron-electron interactions in the 
device are accounted for by the interaction Hamiltonian 
Hi. Its particular form depends on the application; for 
example a common choice in the context of quantum dots 
is given by 



H 



(3) 



Here, the matrix U m i accounts for intra and inter-dot 
interactions, if the labels m and I belong to the same 
or different quantum dots, respectively. The reservoir 
Hamiltonian 



Hr = Y £ ^k(t)bl k b ak . 



(4) 



a£L,R k 



describes non-interacting electrons in the left (a = L) 
and right (a = R) contact. The respective tunneling 
Hamiltonian is given by 



T kl b ak 



(5a) 
(5b) 



with Tjg denoting the couplings between device and reser- 
voir. The operators b^ ak and b a k are electron creation and 
annihilation operators for reservoir states. The abbrevi- 
ations (|5bj) are sums of reservoir operators which will be 
convenient later on. In general, as emphasized by the 
time arguments, all levels and couplings may explicitly 
depend on time through the variation of source, drain or 
gate voltages. 

In order to simplify the notation we will use a com- 
pact multi-index notation in the following. Thereby, the 
multi-index a = (00,01,02) encloses the operator type, 
either annihilation or creation operator, ao = +, — , con- 
tact label a\ = a — L,R, and system state a2 = I. A 
sum over the multi-index J2 a corresponds to a threefold 
sum J2 a Ylf Further, the system operators will be 
denoted by S* a and the reservoir operators, as introduced 
in Eq. (|5b]l. by B a with 



R (+ _ R t 

<? (+) - r, 



B 



(-) 
ad 

(-) 



B a l 
C l ■ 



Therewith, the tunneling Hamiltonian (JSJ becomes 

HsR — ^2 a 0^a'S'a ■ 



(6) 



(7) 



The component ao accounts for the fermionic character 
of the operators. Throughout the article we set ft = 1. 



II. REDUCED DENSITY MATRIX 

The evolution of the full density matrix g(t) is given by 
the Liouville-von Neumann equation, which in the inter- 
action picture with respect to H$ + Hr reads [U, HB| 



[H sn (t),g(t)}_=C(t)g(t) 



(8) 



In the rightmost expression we introduced the super- 
operator £(t) for the commutator with the tunnel- 
ing Hamilton operator. Back transformation to the 
Schrodinger representation is achieved by the super- 
operator 



Ua(t,T)» = U (t,r) • C/t^r) 



(9a) 



with 



U (t,T)=Texp(-iJ^ dt' [H s {t') + H K (t')]\ . (9b) 

Considering only the degrees of freedom of the device, it 
is convenient to define a reduced density matrix a(t) = 
Tre {f?(i)} in the usual way. Starting from Eq. (JSJ) an 
equation of motion for the reduced density matrix a(t) 
might be found in terms of cumulant expansions. Here 
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we will use the partial cumulants, which yields a time- 
nonlocal generalized QME [37l - l40| . The respective equa- 
tion of motion is 



t 

ig- t cr(t)=-i J dr JC(t,T)a(r) 



it, 



OO p 

-i^(-l)"- 1 / dr^ 2n \t,T)a(r) . 



(10) 

The memory kernel IC(t, r) has been expanded in partial 
cumulants with the 2n-th order memory kernel given by 

& 2n \t,T) =fd n ... J dT 2n ^ 2 (11) 
r r 

x (£(t)£(ri)...£(r 2 „_ 2 )£(r))p C 

and the partial cumulants (. . .) pc being defined via a re- 
cursion relation [33] ■ Note, that in Eq. (fT0|) we have an- 
ticipated the vanishing of the odd order cumulants due 
to our choice of the initial density operator, which will 
be discussed below. Obviously, the 2n-th order mem- 
ory kernel contains n super-operators C and therefore n 
times the tunnel coupling i?sR- It has been shown [Il[ 
how this expansion can be connected to a diagrammatic 
description [42}, which provides a clearer interpretation 
in terms of tunneling processes. However, for practical 
calculations one usually considers only the second order 
QME. This is mainly due to the complicated structure of 
nested integrals in Eq. (fTT|) . In the following we will pro- 
vide explicit expressions of the second and fourth-order 
kernels and in next section we will show how an auxiliary 
operator method can be utilized to deal with the integrals. 

A. Second and Fourth Order Memory Kernels 

Firstly, we only consider the second order and fourth- 
order kernels, which explicitly read 

JC^(t,r) = (C{t)C(r)) , (12a) 

* Tl 

/C (4) (i,r)= J d n J dT 2 (£(t)£(n)£(T 2 )C(T)) 

T T 

-(£(*)£(ti))(/:(t 2 )z;(t)). (12b) 

Here the averages (. . .) are meant to be taken with re- 
spect to the equilibrium reservoir state, i.e. (...) = 
Tr r {. . . /9b (to)}- Notice that in writing Eq. (|10[) we have 
assumed that the initial full density matrix factorizes, i. e. 
g(to) — Ps(^o) fJ (io), where the reservoirs are taken to be 
in equilibrium (grand canonical ensemble). Therefore, 
taking the initial time to the infinite past, the QME de- 
scription corresponds to the scheme introduced by Caroli 
and co-workers I34J1 . 



The second order memory kernel in terms of system 
and reservoir operators is explicitly given by the following 
commutator l3lL kSol 



/c< 2 )(v 



E 

a,b 



5 a (t),Cab(i,r)5b(r). 



- • S b (T)C ba (T,t)}_ , 

where we used the reservoir correlation functions 
Cab(*,r) = (B a (*)-Bb(r)) • 



(13) 



(14) 



Due to the fermionic character of the B-operators, there 
are only two possible types of correlation functions 



.(+-) 

ami 

.(-+) 



(t,r) = (Bt m (t)B a i(r)) 
(t,r) = (B al (t)Bi m (r)) 



(15a) 
(15b) 



which we can identify with the lesser and greater tunnel- 
ing self-energies appearing in the NEGF formalism, cf. 
Eqs. (|B4j) . Note that all correlation functions with op- 
erators from different reservoirs vanish since there is no 
direct coupling of the reservoirs. 

For convenience we introduce two correlation super- 
operators, 



C ( sb(t,r). =C ab (t,r)S h (r) 



- • S b (T)C ba (T,t) , 

C ab (t,T)S b (r). 
+ • Sb(r)C ba (T,i) . 



(16a) 



(16b) 



The first definition allows for writing Eq. (|13j) in a com- 
pact form. According to the two correlation functions 
(TTSI) . there are also two types of non- vanishing correla- 
tion super-operators (|16al) . The second definition (I16bp 
is introduced here for completeness. It will be used for 
the fourth-order kernel, which we discuss now. 

The fourth order is more complicated. We summa- 
rize all details in appendix [X] and give here, by means 
of Eqs. (|A6|) and (|A7|) . only the final expression for the 
fourth-order kernel 



K^(*,t)ct(t) = 



,b,c.d 



S a (t), / dri I dr 2 



C atdb(*> r l> T 2> T )< 7 ( T ) 



(17) 



with a new super-operator 



Cdb(^i,r 2 ,r) =A^(n,r 2 )C^(t,r) 



^l(t,T2)C^(n,r), , (18) 



where we have used the correlation super-operators (|16[) . 
It should be noted that each of the sums in (fT7| runs over 
a three-fold multi-index, which comprises the operator 
type, the reservoir and the system state. 
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B. Electric Current 

In order to get a complete description of the nanode- 
vice dynamics an expression for the time-dependent elec- 
tric current is mandatory. In the following we will show 
that such an expression may be found in terms of re- 
duced quantities. The electric current through tunneling 
barrier a is given by the rate of change of the particle 
number in the respective reservoir [2^, [3|| , 

J a {t) = '^(N a ) = -ie{[H,N a }_) 

= -2 e Im j^Tr (i&q e (t))J 

= 2eRe |Tr s Tr R ^ (iBt,(t)q(*)(?(f)) j , 

(19) 



with the number operator for reservoir a given by N a — 

Yl b a ub a k- In the last line we have used the invariance of 

k 

the trace against unitary transformations to switch to the 
interaction picture with respect to Hq. From the equa- 
tion above one sees that the electron current is given by a 
block of the total one-electron density matrix describing 
the coherence between reservoir and system. The deriva- 
tion follows the arguments given for the expansion of 
the memory kernel in terms of partial cumulants, which 
makes use of projection operators It can be shown 
that the current is given by the following expression (43j 



J a {t) = 2e Re I Tr e 



Z Z '/ 2,1-3 

P P P 

^{-lr-'jdr Jdn... j dr 2n ^(H^ a (t)C(T 1 )...C(T 2n ^)C(T)) pc a(r) 



(20) 



r 



with #SR,a(<) = Ei-BL(*) C «(«)- The integrand is very 
similar to the expanded memory kernel (fTTj) . which will 
be used in the next section to calculate the current from 
auxiliary density operators. 



III. AUXILIARY OPERATOR METHOD 

In the following we device an approach which allows for 
the efficient propagation of the QME for the reduced den- 
sity matrix a, cf. Eq. (1101) . We introduce auxiliary den- 
sity operators (ADOs) which contain a time integration 
of the memory kernel operators K.^ and K.^ applied 
on a. Most importantly, these ADOs depend only on 
one time argument which is why they are much easier to 
handle numerically. Their time evolution is given by one 
equation of motion. 



A. Auxiliary density operators (ADOs) 



Using the C super-operators (|16al) we define second 
order ADOs 



n i 2) W=E / drU s (t)C^(t,r)a(r) 

1~ ** 



(21) 



which contain information on the coherence between sys- 
tem and reservoirs as we will show in Sec. IIII El Notice 



that the back transformation lis to the Schrodinger rep- 
resentation is also contained in Eq. ([2"Tj) ; the respective 
super-operator Us is defined as in Eq. © but with Hs 
only. By means of these ADOs, the second order time- 
nonlocal QME in Schrodinger representation can be com- 
pactly written as 



i^a(t) = [H s ,a(t)]_-iY l [Sa,^(t) 



(22) 



This form is identical to the standard one [3l|, |36( , which 
becomes clear if one of the ADOs is written explicitly, 
e. g. a = (+, a, m) 



t 

U^(t)=J2jdr 



(23) 



(c£-\t,T)U S (t,T)4o-(T)Ul(t,T) 

-CLVi^t) Us(t^)a( T )4ul(t,T) 



Before discussing the propagation of the equations for the 

reduced density matrix o~(t) and the operators Ila (t) we 
will derive the equations for the next higher order. 

In analogy to the second order case we introduce a 
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fourth-order ADO 



usual convention, i.e. 



I I Tl 

ni 4) W = J2 J dT J dTl J dT ^ u s(t) 

b > C > d t T T 

x [^c(T 1 ),ci 4 c ) db (i,T 1 ,r 2 ,r)a(7 



(24) 



which allows for a compact form of the fourth order kernel 
([T7| . Thus we can write the QME up to fourth order in 
Schrodinger representation as 



i|a(i)= [H s ,a(t)l 



(25) 



with the Il( n ) operators defined by Eqs. (gU) and (12H) . 
Having this compact form we may now derive equations 
of motion for the ADOs in the so-called wide-band limit. 



B. Wide-band Limit (WBL) 

The proposed propagation method is based on the ob- 
servation that a closed set of equations of motion for 
a and IT 2 ) can be obtained, if the reservoir correlation 
functions (|15p can be expanded in a sum of exponentials 
[27l |44| . In particular, this is possible if the so-called 
level- width function, 



(26) 



can be approximated by a sum of Lorentzians [44| . If T a 
is only slowly changing with e on the energy scale given 
by the system, one can employ the wide-band approxi- 
mation, i. e. assuming an energy independent level- width 
function, T a (e) = T a . For the rest of this article we 
will concentrate on this limit. Using an expansion of the 
Fermi-Dirac function as described below, one still finds 
the correlation functions to be sums of exponentials plus 
one term being proportional to a <5-function in the time 
domain, which is a direct consequence of the WBL. To see 
this explicitly, we consider the calculation of the reservoir 
correlation functions (|14[) . which involves an integration 
over the Fermi function, see Eqs. (1151) . 



C ab (t,T) = (c +iffRt B a c- iffRt c +iffRT B b c- iffRT ) 

The sign in the Fermi function / and in the exponen- 
tial characterizes the type of the correlation function, 
namely ao = + refers to C' H ^ and ao = — to ^, cf. 
Eqs. dT5|) . As usual /3 denotes the inverse temperature. 
The level-width function as defined in Eq. (j2l))) also de- 
pends on the order of the bath operators. We follow the 



ab 



T a for a = — 
r£ for a = + 



(27) 



For convenience we also neglect an explicit time- 
dependence of the tunneling matrix elements. Then, in 
the WBL the correlation functions read 



C ab (*,T) =T° 



ds 
2^ 



/(a /3(e-^))e aoi£(t - T) • (28) 



In order to perform the energy integration we expand the 
Fermi function in terms of a finite sum over simple poles 



1 



1 + e ± / 3 ( £ -' J ) 
1 l"> ' 



(29) 



with xtp = + x p //3 = (xap)* ■ Instead of using the 
Matsubara expansion (45|, with poles x p — i7r(2p— 1), we 
use a partial fraction decomposition of the Fermi func- 
tion [331 ] , which converges much faster than the standard 
Matsubara expansion. Furthermore it allows to estimate 
the error made by truncating the sum ([29]) at Np terms. 
For this decomposition the poles x p = ±2y^ are given 
by the eigenvalues z p of an Np x N-p matrix |33j | . We take 
the branch of the root yjzp~ such that Im (x p ) > for 
all p. Thus all poles Xp (Xp) are in the upper (lower) 
complex plane. 

Finally, one obtains by means of Jordan's Lemma the 
following expressions, 



C ab (i, r) = -r Bb <f(*-r) - J2 C *b P (t, t) , 

P =i 



C ahp (t, t) = — TahC 1 

C a b P (i, t =F 0) = "^r a b 



(30a) 
(30b) 

(30c) 



with s = sgn(i— t) and Xnp = a oX^p- Apart from the 
first term in Eq. (|30ap the correlation functions are sums 
over exponentials. The same strategy might be used 
if the level- width function T a (e) is given by a sum of 
Lorentzians. 



In this case the first term in Eq. pOal) is 
also a sum of exponentials; see Ref. [13, [lH, |44| . 

If the energies in the reservoir Hamiltonian ^ arc 
time dependent the expression in the exponential in Eq. 
([28]) has to be replaced by an integral, aoie(t — r) — > 
/' dt'aoi(e + A a (t')). Consequently, the auxiliary modes 
take the form, 

xZ P (t) = [Ma + A«(*)]+ Xp/P- 

In the following, we will suppress the time-label to facil- 
itate readability of the equations. 
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The expansion (|30|) allows us to find the equations of 
motion for the reduced density matrix and the second 
order ADO. Moreover, we will show that one can also 
obtain an equation of motion for the fourth-order ADO, 
which eventually gives corrections to the second order 
equations, but also describes new processes. 



C. Equations of motion for the second order ADOs 

Using the expansion of the correlation functions given 
by Eqs. (|30l) we may rewrite the correlation super- 
operators (Unj) as follows, 



1, 



4&(t,T)» =-T ab [S b (T),.],6(t-T) 



(31a) 



(31b) 



where the sum runs over the auxiliary modes. The partial 
super-operators have a very simple structure, 



Cg,(*.r)' 

1(2) 



C ( Z(t,r) [S h (r), 

(2) 



A K :L{t,T). =C^(t,r)[S h {r), 



(32a) 
(32b) 



abp V ) / abp \ 

which follows from Eq. (|30bl) . 

Using the expanded super-operators (|3T|) in Eq. (f2l"j) 
readily leads to an expansion of the second order ADOs 
in terms of partial operators, 

= 2E r ^b,a(t)L -XXpW • (33) 



(2) 

The partial operators H ap (t) are similar to the definition 
(f2Tj) but with replaced by . Their equations 

(2) 

of motion define the time evolution of H a (t) . These 
equations read 



d_ 

~d~t 



E4 r ab [S h ,a(t)}_ 

b p 



ng?(*)> 



(34) 



whereby each of them contains one auxiliary mode x ap . 
The initial condition can be found from the definition of 

( 2) 

the partial operators, H ap (to) = 0. 



D. Equations of motion for the fourth-order ADOs 

Since the fourth-order ADOs [Eq. (|2"4"|) ] are given in 
terms of the super-operator we first need an expan- 
sion for these objects. To this end we apply Eqs. (|3"Tj) 



in two steps. First we expand only correlation super- 
operators containing time label t. This gives 



C 



(4) 

acdb 



(t,T 1 ,T 2 ,r)> 



4d(n,r 2 )-r ab [S h (r),.]_S(t-r) 



E C acdbp(^ T l' T 2^)» 



5(t-T 2 ) 



(35) 



with 



Cilpt^l, T 2 ,T) = A%(T U T 2 )C£ p (t,T) 

j(2) 
adp 



,(2) 



^ht^C^^r). (36) 



With these results one can expand the fourth-order ADO 
given by Eq. Notice that the nested time integrals 

in are such, that the terms in Eq. (pi5j) containing a 
(^-function vanish. Therefore, one has 

ni 4 )(t)= E n $(*) ( 37a ) 
p 

* t Tl 

ngK*) = - E j dT J dTl J dr 2^s(t)x 



b.c.d 



to 



[-5 C (n), 



■(4) 

acdbp 



(t,Ti,T 2) r)cr(r) 



(37b) 



In analogy to the second order case one can give an equa- 



tion of motion for Tl a 4 J , 



a 



[ff s (t),n$(t)]_ -xa P ni 4 p )(t), 



(38) 



which contains only one auxiliary mode. The last two 
terms result from the time-derivative of the system prop- 
agator and the reservoir correlation function, respec- 
tively. The first term, denoted by $>*c P (t), comes from 
the derivative of the integrals and hence is given by 



b,d r _ 



Us(t)C^ dbp (t,t,T 2l T)a( T ) 



(39) 



The equation of motion is supplemented by the initial 
condition found from Eq. (|37bl) to be U^(t ) = 0. Now 
we expand the remaining two correlation super-operators 
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in C' 4 ) formerly involving time label t\, 

C^bp(*,*,7a,T)« 



- 2 Fcd 

,(2) 



S d (t),cW p (t, T ).] + 5(t-T 2 ) 
1, 



-E C icdb P9 (*'*'^,r). . (40) 

Inserting this expression into $icp(i), one observes that 
due to the integral boundaries the term containing the 
delta function S(t — r) vanishes and one is left with 



^) = 7EM^ d ' n ^ (<) j -E $ acUM> (41) 
d + q 

and 



z z 

J dT J dT 2 UsW 



(4) 

acdhpq 



(t,t,T 2 ,T)a(T) 



(42) 



Obviously one gets two qualitatively different contribu- 
tions to the equation of motion for II^ 4 ) . The first one be- 
ing proportional to IT^ 2 ' contains only one excitation (ad- 
ditional electron or hole) in the reservoirs and gives cor- 
rections to the equation of motion for II^ 2 ). The second 
contribution, given by $itp 9 , describes two excitations in 
the reservoirs, i. e. the onset of genuine co-tunneling. 



In summary, we write for the equation of motion of the 
fourth-order partial operators, 



i|ni 4 ) W =-i£ 



Sc,-E r cd \s dl n$(t)} -E^tLW 



fls(*),n<i>(*) 



(43) 



In order to get a complete description we also need the time evolution of $icpq, which can be obtained by using again 
a differential equation, 

ij^itU*) = 1 E (4dU*' * + ^^P W - 42p(*> t + 0) n $(*)) + $ it, (*)] _ - (Xap + Xc 9 ) $W,(t) ■ (44) 



The initial condition is easily obtained, $itpq(^o) = 0. 
Using Eqs. (1321) and (1421) one can also show that $itpg is 
a traceless operator, 



Tr S {$acp 9 (t)} =0 



(45) 



This property will be important in the next section, when 
we consider expectation values for example to calculate 
the electric current. 



E. Electric Current 

So far we have given a closed set of differential equa- 
tions for the reduced density operator and newly intro- 
duced partial operators. By means of these quantities 
one can calculate the time evolution of any observable of 
the system. Moreover, one can also compute quantities 
like the time-dependent electric current, which do not 
belong to the system state space alone. 



Up to fourth order the electric current (f2"0")) becomes, 
jW(i) = 2eRe j^Tr s [c m n( 2 m +)(t)" 



E^< 



c n (4 ' +) 



(0 



(46) 



where the definition of the ADOs given by Eqs. (|21[) and 
(PHI) has been used. Using the expanded forms of the sec- 
ond and fourth order operators, cf. Eqs. (|3"3"f and (|37al) . 
one finds 



Tr s [ Cm ni 2 / +) (t)- Cm nit +) W 

= Ez^ ( * )_cr ^' ( * ))r 



cc,jl 



E Tr skni1 P +) W + ^n^(t) 



(47) 
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where <J m j{t) — Trs{cjc m er(i)} is the one-particle and 

&mj(t) = Trs{c m cjer(i)} is the one-hole density matrix 
of the system. The current is therefore given in terms 
of the reduced one-particle density matrix and the aux- 
iliary operators. For the last terms involving the sum 
over expanded ADOs one might now use the equation of 
motions derived in the previous section, cf. Eqs. (|34p and 
©. 

In an occupation number representation the traces in 
the current expression as well as the equations of mo- 
tion can be solved numerically for an arbitrary device 
Hamiltonian ([2]) including interactions. As an example 
we will show results for an interacting double quantum 
dot system in Sect. [V] 



IV. INFINITE ORDERS 

In the previous section we have considered the cumulant 
expansion up to fourth order. Naturally, the question 
arises how the higher orders influence the dynamics. In 
principle, one has to sum up all orders to get an exact 
description. This problem has been tackled, for example, 
diagrammatically [42j or by using a path-integral formal- 
ism [28| . Moreover, the connection of cumlant expansions 
obtained from the projection operator technique to the 
diagrammatic expansion has been discussed recently [4l| . 
Other approaches consequently use many-body states to 
derive equations of motion [2g, H|| [47| ■ 

In the following, we will first show for non-interacting 
electrons that summing up all orders with one excitation 
yields a QME, which resembles the equations of motion 
from the well-known NEGF formalism, cf. appendix [5] 
It contains level broadening of the system states due to 
coupling to the reservoirs. Further, we adopt the path- 
integral formalism [28| and device a propagation scheme 
for systems with electron-electron interactions similar to 
the one developed in the previous section for the WBL. 
We restrict the discussion to single-electron excitations 
in the reservoirs only. 



A. Noninteracting Electrons 

Neglecting electron-electron interactions we can com- 
pare the expressions obtained from the QME and the 
NEGF formalism. For example, the electric current up 
to fourth order is given by 

j£\t) = 2e Re £ J i £ {a ml (t) - a rnl (t)) T a , lm 

- V ( p {2 ' +) (t) + p^+i (t)) 1 

P ) 

(48) 



and therefore in terms of single-particle quantities a m i (t) 



andP^ ±) W = Tr, 



with n = 1, 2. In 



order to get a complete description one has to calculate 
the equations of motion for these matrices, which may be 
obtained directly from the respective QME. 

Considering the matrices P^^t^ we find from Eq. (|34|) 
that their time derivative for the second-order case is 
given by 

d_ 

-w 



4^ (*) = - \ E t^w + r 



Yh mi p {2 . 



h P^^'(t) - v+ P 



ajl 
(2, 



ap ap\ml 



(*) 



(49) 

with h m j = e m S m j + V m j as an abbreviation for the de- 
vice Hamiltonian (j2J|. The fourth-order contribution is 
obtained from Eq. (|4"3"]l and reads 



i A p( 4 <+) a) - - \- V r / p {2 ^> (t) 

1 Q^ r ap;ml\ L ) ~ L n / , 1 Q ' ™3 r av.il \ b ) 
a'j 



mj^ap-Jl y L ) 
(4, + -) 



X, 



+ 



>(4,+: 



ap ap;ml 



al <x'm,pq 



(0 



(t) 

(50) 



whereby the term Trs {^it^a'z pg(*)} ' s easu y seen to 
vanish by using Eqs. f4"2")l and (f^Tj) . 

Comparing Eq. (|53|) with the exact NEGF result given 
by Eqs. (|B1[) and (IB6p suggests that our results constitute 
the first terms of an expansion eventually leading to 



P 



(+) 



cx.p\ral 



(*) 



Ep(2n,+) 
ap;rnl 

n=l 



a ' 



(51) 



where each term p( 2 (™ +1 )'+) gives rise to corrections for 
the lower-order term p( 2n >+). Indeed, one sees that gen- 
eralizing Eq. (f5T))) to 

. d 



_ >(. 2n >+){+\ - - -X^V p(2n-2. + ), 

L ap;ml 2 < m 3 awM 1 

a'j 

+ s Th -p^+h 

j 

,(2n,H 
^ap* ap;mi 



Aop J av.mi \ / 



(52) 

with n > 2 yields the correct equation of motion for 
Pap-mi- By comparing with Eq. (|48|) the current is given 
by ' 



J a (t) 



2e Re ^ \ J ^ (a ml (t) - a ml (t)) T c 

m L I 



E 

p 



p (+) ft) 

ap-mm \ J 



(53) 



This expression is identical to the respective result ob- 
tained from the NEGF formalism, cf . appendix [B] 
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B. Interacting electrons 



noted by auxiliary mode p. Its equation of motion reads 



Inspecting Eq. (j43]l . which led to Eq. ([50| in the first 
place, one is inclined to conclude that for arbitrary order 
n > 2 one has to use a modified EOM, 



dt ap {) 



So, 



t£[Sc,*^(t)] ■ ( 54 ) 



This equation generalizes Eq. (|43| and leads to the cor- 
rect form given by Eq. (|52[) . Therefore, this generalized 
EOM yields the exact dynamics for vanishing electron- 
electron interactions. 

For a rigorous derivation one has to sum all higher- 
order cumulants. Instead of explicitly summing all 
higher-order cumulants in Eqs. (fTU|) and (I2U1) we use the 
hierarchy of equations of motion for the ADOs derived 
elsewhere 28]. Since we are interested in the WBL we 
have to modify the equations given therein. In particu- 
lar the term containing the ^-function in the expansion 
for the correlation functions given in Eqs. (|30l) leads to 
a substantial simplification. The equation of motion for 
the reduced density matrix is almost identical to Eq. ([25]) , 



i^a(t) = [H s (t),*(t)] 



iJ2 [S«,fia(*) 



(55) 



The new operators II a involve all orders Il( 2n ) . Using the 
expansion of the correlation functions given in Eqs. (|30|) 
we obtain the following expansion 



na(i)=X>a n) W 

n=l 

= Ei r *b[Sb,<7WL-Ena P W • (56) 



i|n ap (t) = 



[H s (t),U ap (t) 



XapIIap(^) 

Sd, n ap (<) 



1 y t h^c! n acpg (i) 

cq 



(58) 



where we recognize the structure of Eqs. (|3"4")l and (|4"3"j) . 
Obviously, Eq. (j58|) contains couplings to the Oth and 
2nd level of the hierarchy, which are given by the first 
and the last term, respectively. The operator n acpg = 

~Y^=i ^acpq surmounts all contributions with two aux- 
iliary modes. The second level of the hierarchy is then 



given by the EOM for II 
hierarchy, we set II 



acpq 



In order to truncate the 



0, which does not influence 



acpq • 



the results for non- interacting electrons (see Eq. ([50)) ). 
This Ansatz is not unique, however, it is very convenient 
since it preserves the linearity of the EOMs. The differ- 
ential equations (f5"5"j) . (|5H|) and the expansion (151))) pro- 
vide a closed description of the evolution of the reduced 
density operator. These equations will be denoted as ef- 
fective quantum master equations in the following. As it 
was shown above, these equations are sufficient to calcu- 
late single-particle quantities like the electron current for 
vanishing electron-electron interactions. This is true for 
non-interacting electrons, but also for a very large inter- 
action strength, where only one electron can enter in the 
system. In contrast to the finite order equations (second 
and fourth order QMEs) , the effective QME consistently 
treats the broadening of the electronic levels due to back- 
action of the electron reservoirs on the system, which is 
also included in the NEGF formalism. The correct han- 
dling of the level-broadening is a necessary ingredient for 
the description of electron transport as we will show in 
the next section for a specific example. 



NUMERICAL RESULTS: DOUBLE 
QUANTUM DOT 



which reminds of Eq. (|3"3")l . The electric current is then 
calculated from the single-electron ADOs only [28| 



J a {t) = 2e Re ^Tr< 



»n+,(t) 



(57) 



which resembles Eq. (|20j) . 

The first level of the hierarchy is given by the EOM for 
the auxiliary operator n ap , which describes the single- 
electron partial ADO with one reservoir excitation de- 



In this section we will implement the propagation 
schemes developed in the previous sections. In particu- 
lar, we will investigate the influence of higher-order terms 
in the QME and compare the results with other meth- 
ods. To this end we consider a double quantum dot 
(DQD) system, where the two dots are coupled in se- 
ries and each dot is connected to an electron reservoir. 
Over the past years this system has become a paradigm 
for electron transport through a nontrivial quantum sys- 
tem [3(3, Ei, 48]. Moreover, the advent of explicitly time- 
resolved experiments, such as coherent control investiga- 
tions allows to test the suggested theories. 
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A. Setup 

In following we will consider only spinless electrons, 
which simplifies the calculations to some extend, but is 
not mandatory otherwise. The Hamiltonian of the reser- 
voirs and the coupling of the DQD to the reservoirs is 
given in Eqs. ((!]) and ([5]). The Hamiltonian describing 
the DQD is explicitly given by 



^bias I 



K- IT] 

bias L J 



+ V(c\c r + c\ci) + Uclciclcr , 



(59) 



where the localized states of the left and right dot are de- 
noted by a £" and "r" , respectively. The first three terms 
in Eq. ([59)) describe an effective two-level system, which 
is characterized by the interdot tunnel coupling V and 
energies +e and —e, respectively. The interdot Coulomb 
repulsion-strength is denoted by U. The eigenenergies of 
the Hamiltonian (l59l) are 



E = 0, Ef = ±y / e 2 +V 2 , E 2 = U. 



(60) 



These energies correspond to configurations with zero 
{Eq), one (Ef) and two (E 2 ) electrons in the DQD, re- 
spectively 

Since the two dots are coupled in series the level- width 
functions in the WBL contain one non-zero element, 



l — 



r/2 o 





R = 





o r/2 



Here, we have assumed a symmetric coupling of the left 
and right dot to their respective reservoir. Following 
Ref. [26| we assume a symmetric voltage drop across the 
system, /i L = -/x R = Vbias/2 und e = r/2. 

The QMEs are conveniently represented using the 
occupation number basis consisting of four states: 
|0),|^),|r) and |2). Consequently, all operators become 
4x4 matrices. The matrix equations resulting from the 
equations of motion are prop agated using a fourth or- 
der Runge-Kutta scheme [49(. The number of auxiliary 
modes is = 120. In all cases we start the simula- 
tion by suddenly coupling the initially empty DQD to 
the reservoirs at t = 0. The temperature of the electrons 
in the reservoirs is set to T = 0. IT/ fee- 

We consider the second and fourth order QMEs, here- 
after referred to as QME2 and QME4, given by Eqs. ([25]). 
(IM| and (g3D, where for the latter we set $$ tPq = 0. The 
QME resulting from the hierarchy will be denoted as ef- 
fective QME (QMEe). We will also present results for 
the NEGF method and the Markovian QME, which are 
explained in Add. [Bl andlCl respectively. 



B. Stationary Current 

Before turning to the transient behavior of the differ- 
ent QMEs we will address the stationary current as a 
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FIG. 1: Stationary current Jl vs bias voltage Vbias for differ- 
ent values of the Coulomb repulsion U. Left panel: results 
obtained from the effective and the Markovian QME. Right 
panel: results obtained using the second order and fourth or- 
der QME and the effective QME. Vertical lines indicate the 
transition energies. Arrows show the values obtained inde- 
pendently in the large- bias limit [50l. l5l|. 



function of the bias voltage. To this end we propagate 
all QMEs as described in the previous section up to a 
final time t = 60T' 1 . 

Figure [1] shows the stationary current for different val- 
ues of U. The left column provides a comparison of the 
effective QME with the Markovian QME. For U = we 
also show exact results using the NEGF formalism [35| 



and the large-bias limit [50j, [51( . The right column gives 



a comparison of QMEe with the second and fourth order 
counterparts. 

The current displays a single (U = 0, 16r) or multiple 
step behavior (U = 4r), which is monotonic for the ef- 
fective and Markovian QME, but non-monotonic for the 
second and fourth order order calculations. The position 
of the steps is best understood in terms of the eigenener- 
gies (IM|) of the system Hamiltonian and the allowed tran- 
sitions between them. Transitions are only allowed by 
changing the electron number by one, i.e. either Ofsl^ 
or 1 ± o 2. The steps in Fig.Q]occur, when a new transi- 
tion energy, Eab — Ea — Eb, becomes accessible within 



11 



the transport window {/j,l,Vr} = {^bias/2, - V hia , s /2}. 
Therefore, a step in the current- voltage curve is expected 
for Vbias/2 = ±Eab- The transition energies are given in 
Tab. U and the respective voltages are indicated in Fig.Q] 
by vertical lines. Due to the symmetry of the setup only 
up to three transitions are visible in Fig. [TJ 

The Markovian results show clear steps for all values of 
U. The width of the rising edges is of the order UbT. In 
contrast to that, the effective QME yields much broader 
steps. This is a result of the broadening of the energy 
levels due to the coupling to the reservoirs, which is not 
present in the second order description underlying the 
Markovian approximation. In the present case the cou- 
pling strength to the reservoirs is larger than the tem- 
perature, r/2 > k^T, and thus the rising of the step is 
mainly determined by the tunnel coupling. For a suffi- 
ciently large bias voltage all energy levels are within the 
transport window regardless of being broadened or not. 
In this limit both methods yield the same results, i.e. the 
Markovian approximation is sufficient. This is readily 
seen in Fig. [1] The small deviations for finite U might 
be attributed to neglecting the contributions of the two- 
electron ADOs. 

The second and fourth order calculations show promi- 
nent features at the transition energies, which are not 
expected within the usual picture of sharp energy lev- 
els. Figure [2] shows the deviations of the QME2 and 
QME4 results from the effective QME calculations for 
U — AT. The behavior close to the transition energy is 
similar for both QMEs. This behavior may be traced 
back to an energy renormalization which is determined 
by the imaginary part of the unilaterally Fourier trans- 
formed correlation functions [46|, [52j , 

oo 

CabM = J dt'C ab (t')e X p[±iu>t'} . (61) 



In Fig. [5] the behavior of ' is exemplarily shown as 
a function of the bias voltage. One observes pronounced 
features at the transition energies. In the context of non- 
Markovian dynamics this can be related to the Lamb- 
Shift Hamiltonian, which may lead to deviations from 
the canonical stationary state (53[. Explicit time propa- 
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TABLE I: Transition energies for allowed transitions. The 
energies are calculated by using Eqs. (|60[) with e = r/2 und 
V = T. 




^bias ^ 1 

FIG. 2: Upper graph: Laplace transform of the reservoir cor- 
relation function according to Eq. (|61[) . Lower graph: 
Stationary current Jj, vs bias voltage Vt>i as for Coulomb repul- 
sion U — 4F. Curves show the difference with respect to the 
effective QME results. Vertical lines indicate the transition 
energies. Deviations for the Markovian QME are indicated 
by the gray-shaded areas. 

gation schemes intrinsically account for the renormaliza- 
tion, which has been discussed in the context of bosonic 
reservoirs [32| • The non-monotonic behavior observed in 
the present case indicates the breakdown of the cumulant 
expansion with a finite number of terms, which is seen 
to be valid for couplings smaller than the temperature 

r<teT0, 

C. Sudden Switching 

The strength of the formalism given in Sec. IIIIBI lies 
in the ability to study transient responses of the sys- 
tem. Here we consider the time-resolved occupation and 
current after suddenly coupling the initially empty DQD 
to the reservoirs. Figure [3] shows the results obtained 
from the different QMEs for non-interacting electrons at 
^bias = 3r. In all cases one observes an initial tran- 
sient response to the sudden coupling and the eventual 
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t[l/T] I HIT] 

FIG. 3: Time-resolved occupations n^, r and currents Jl,r for 
U = 0. At t = the initially empty DQD is suddenly coupled 
to the reservoirs. Left (right) panel shows occupation of the 
left (right) dot and the current through the left (right) barrier, 
respectively. The bias voltage is Vbias = 3I\ 

attainment of a stationary state, which was discussed in 
the previous section. Comparing the QME results to the 
exact NEGF calculations one observes that the effective 
QME indeed yields the correct dynamics in accordance 
with Sec. IIV Al In contrast, the second and fourth order 
calculations reveal substantial deviations from the exact 
results. In particular, we observe a slower damping of 
the oscillations, which is expected since the additional 
term in the effective QME provides a damping for the 
ADO. This term is not present in the second oder case 
and only partially accounted for in the fourth order case. 
The deviations in the long-time limit correspond to the 
discussion for the stationary case. 

Going to larger bias voltages one finds that these de- 
viations vanish and eventually the Markovian result is 
obtained. This is illustrated in Fig. [4] The difference be- 
tween the Markovian approximation and the other meth- 
ods for t — > are due to the lack of memory in the Marko- 
vian approach. This behavior complements the discus- 
sion of the large bias limit in the previous section. 



VI. SUMMARY 

We have developed a propagation scheme for the re- 
duced density matrix by means of a time-nonlocal QME. 
Hereby, we obtained expressions for the memory kernel 
up to the fourth order, which has not been discussed 
before. In order to get a viable numerical scheme we in- 
troduced auxiliary density operators. An expansion of 
these operators in terms of auxiliary modes, given by a 
partial fraction decomposition of the Fermi function [33| , 



yields a set of coupled differential equations that can be 
solved using standard methods. To investigate the valid- 
ity of the fourth-order equations we considered the par- 
ticular case of noninteracting electrons. A comparison 
with exact results for the electric current obtained from 
an NEGF formalism [22| showed that the fourth-order 
description is not complete. Based on this analysis we 
proposed equations of motion for higher-order auxiliary 
density operators, which leads to the correct dynamics 
of single-particle observables. Interestingly, these equa- 
tions only involve quantities with one auxiliary mode, 
which is interpreted as a single excitation in the reser- 
voirs. The consistent treatment of this excitation up to 
all orders in the tunnel coupling with the contacts leads 
to a broadening of the energy levels, which is essential 
for a complete description of the dynamics. Further, we 
compared the proposed effective QME with the results of 
a path-integral approach [28[ and found that the QME 
corresponds to the first level of the hierarchy. 

Finally, we showed and discussed the results of the 
propagation scheme for the specific example of a dou- 
ble quantum dot. Hereby, we verified that the effective 
QME with a single auxiliary mode is in agreement with 
NEGF results for noninteracting electrons. For a finite 
interaction strength we observed the appearance of un- 
expected structures in the current-voltage curve for the 
second and fourth-order calculations. These findings in- 
dicate the breakdown of the finite-order expansions for 
strong tunnel coupling. 

In total, the combination of auxiliary density opera- 
tors and the auxiliary mode expansion yields an efficient 
method which allows for a numerical description of time- 
resolved electron transport. As we have shown, the in- 
clusion of higher-order terms is a necessary ingredient for 
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problems involving strong tunnel couplings or low tem- and to investigate realistic schemes for coherent manip- 

peratures. This opens the possibility to study the re- ulation of nano-devices. 
sponse of complex systems to time-dependent drivings 

I 



Appendix A: Fourth Order Partial Cumulant 

The partial cumulants may be defined in terms of a recursion relation (37l . l54j , 

(C(t) . . . C(r n ) . . . £(r)) pc = (-1)^ {£(t) • • ■){C{r n ) ■■■){■■■ C(r)) , (Al) 

where the sum is over all partitions and g is the number of partitions in the term. The partitions are such that the 
time-ordering t > . . . > r n > . . . > r is preserved. Taking (£(£)) = the fourth order partial cumulant is 

(£(t)£(r 1 )£(r 2 )£(r)) pc = (£(t)£(n)£(r 2 )£(r)) - (£(i)£(n))<£(r 2 )£(r)) , (A2) 

where the two-point correlation operators are given by (|13p . The four-point correlation operators can be expressed in 
system and bath operators, by expanding the commutators and taking the partial trace, 

(£(t)£(r 1 )£(r 2 )£(r)) • - [s(t), S(n)F . —f . 5(n)] _ , (A3) 

with 

T* = (B(t)B( n )B(r 2 )B(r)) S(t 2 )S(t) . -(B( T )B(t)B( n )B(T 2 )) S(r 2 ) . S(r) 

- (B(r 2 )B(t)B( n )B(T)) S(t) . 5(r 2 ) + (B(r)B(v,)B(t)B(ri)) • S(t)S(t 2 ) (A4a) 
f* = {B( Tl )B(t)B{T 2 )B{ T )) S{t 2 )S(t) . -{B( T )B( Tl )B{t)B(T 2 )) S(t 2 ) • S(r) 

- (B(T 2 )B( n )B(t)B( T )) S(t) • S(t 2 ) + (B{T)B{T 2 )B( Tl )B{t)) . S(t)S(t 2 ) (A4b) 

For self-adjoint reservoir operators B in the expression above one finds f = J 7 ^ (see also [55j . where a derivation was 
given for bosonic operators). 

For the case of non-interacting electrons in the reservoir with an initial thermal state, we can further simplify the 
four-point correlation functions using the fermion version of Wick's theorem [56l appendix 3], 

(B(t)B(n)B(r 2 )B(r)) = {B(t)B(n)){B( n )B(r)) - (B(t)B(r 2 )){B( n )B(r)) + (B(t)B(r)) (B( ri )B(r 2 )) . (A5) 

Using the correlation super-operators ([T5| and Eqs. (|A4[) and (|A5[) one gets 

(£(t)£(n)£(T 2 )£(T)) • - [S(t),C(t,n)[S(T 2 ) i C(T 2 ,T).]_}_ 

+ [S(t),[S(T 1 ),(A(n,T 2 )C(t,T)-A(t,T 2 )C(n,T))»} + ]_ . (A6) 

Moreover, with the definition of the correlation super-operators, (fl~3"|) and (IT^l) . one finds 

(£(t)£(r 1 ))(£(T 2 )£(r)). = [S(t),C(t,n) [S(t 2 ),C(t 2 ,t)»}_]_ . (A7) 



Therefore the fourth order partial cumulant is given by (IA6|) without the first line, which is canceled by the expression 
above. This result was used in Eq. (|17p for the fourth order memory kernel. 

I 



Appendix B: Non-Equilibrium Green Functions approach of Sect. ITVl For a more detailed derivation see 

(NEGF) Ref. [12]. 

The current J a through the barrier connecting lead a 
The description of electron transport by means of non- anc j the system is given by 
equilibrium Green functions has been widely used [57j 
since its first formulation (35j | . Here, we summarize the 
main equations in a form, which differs slightly from the 

usual one, in order to facilitate comparison with the QME J a (t) = 2e Re Tr {II Q (t)} . (Bl) 
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with current matrices 
t 

n, 



a(t)= J rfti [G>(*,* 1 )S<(t 1 ,t) 



-G<(t,ti)E>(t 1 ,t)] , (B2) 

describing the flow of electrons from the reservoir into 
the system and vice versa. 

By means of the Keldysh equation [57j for G > and 
Dyson- type equations [57[ for G a and G r , respectively, 
one obtains as equation of motion for the density matrix 
a(t) = -iG< (t,t) 



iJ^W = [h(t),<r(t)]_ 



i^{n Q (t) + nl(t)} 



(B3) 



which resembles the structure of a QME. 

In order to obtain the time evolution of the current 
matrices H a one can write the self energies 



S<(t 1 ,t) = +iC(+-)(t,ti) 



(B4a) 
(B4b) 



in full analogy to the correlation functions (|30|) as a finite 
sum 

£<->(ti,t)= ±ir a (ti,t)J(t-ii) 



T C a b, p (ti,t) 
P =i 



(B5) 



with S a p(ti,t) = iC ap (ti,i). Using these expansions in 
the definition of the current matrices, cf. Eq. (IB1|) . one 
obtains 

n a (t) = -r a (t 1 t)--tr{t)r a (jt,t)-^2u ap {t) , (B6) 



which defines auxiliary current matrices Tl ap . Using Eqs. 
(|3U1) the equations of motion for H ap read 



i^-U ap {t) =S ap (t,t) 



(B7) 



h(t)-x^i-rr(t,t))n^(t) 



Appendix C: Markov approximation 

Since there tends to be some confusion about the pre- 
cise meaning of the Markov approximation we give in 
the following the definitions used in the present work. 
Starting with the second order QME, cf. Eqs. ([21"]) and 
(|22|) . we first substitute <t(t) by U&(t - r)a(t). This ren- 
ders the QME local in time and the same equation could 
have been obtained directly from the Tokuyama-Mori ap- 
proach (58|. Next we set the initial time t = — oo and 
also set t — t = t' in the integral. The resulting QME 
reads 



ij^(i) = [H s ,a(t)}_-iTZa(t) 



(CI) 



with the usual Redfield super-operator (e.g. see |36j |) 



DC- 

U * = Yj I dt ' [ 5 a,C7 ab (i')SbH')' 



ab 



- . S b (-t')C ba (-i')]_ (C2) 



and Sh(—t') = Us(t')Sb- Further, using the eigenbasis of 
Hs one has to evaluate the Laplace transform C^^Eab) 
of the reservoir correlation functions (cf. Eq. (|o"T1) L which 
can be done analytically for the WBL. The final step 
yielding the Markov approximation consists in neglecting 
the imaginary part of the integral. 
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